function y = cholpsd1(x)
% CHOLPSD1(X)
%   computes a Cholesky decomposition for a positive semi-definite matrix
%   uses CHOLINC

% February 28, 2002
% Sungbae An

% select method
mchol=1;

% the following two procedure make the same output
if mchol
	z = sparse(x);
	y = full(cholinc(z,eps))';
else
	y = x;
	n = size(y,1);
	for k=1:n
		if y(k,k) > 0;
			y(k,k) = sqrt(y(k,k));
		    if k < n;
			    y(k+1:n,k) = y(k+1:n,k)/y(k,k);
				y(k,k+1:n) = zeros(1,n-k);
	            for j=(k+1):n
		            y(j:n,j) = y(j:n,j) - y(j:n,k)*y(j,k);
			    end;
	        end;
		end;
	end;
end